function [HHstates,wage_moments,transitions]=data_simulation_A12_end(Par1,ff_1,wf_1,fi_1,wi_1,ff_2,wf_2,fi_2,wi_2,Wff,Wfi,Wfn,Wif,Wnf,Wii,Win,Wni,Wnn)


df_1=Par1(1);
di_1=Par1(2);
lnf_1=Par1(3);
lni_1=Par1(4);
lfi_1=Par1(5);
lif_1=Par1(6);
df_2=Par1(7);
di_2=Par1(8);
lnf_2=Par1(9);
lni_2=Par1(10);
lfi_2=Par1(11);
lif_2=Par1(12);
p_1=Par1(13);
q_1=Par1(14);
p_2=Par1(15);
q_2=Par1(16);
alphaf=Par1(26);
betaf=Par1(27);
alphai=Par1(28);
betai=Par1(29);
% Data simulation (MSM)

N_households=10000; % # households
maxtime=400;   % # quarters
rng('default'); % set seed

% Pre-allocation
value=zeros(N_households,maxtime);
state=zeros(N_households,maxtime);
wage_draw1=ones(N_households,maxtime);
wage_draw2=ones(N_households,maxtime);

for iter=1:N_households
    
    time=1; % households are born in the situation 9: "NN"
    scale = lni_1 + lni_2 + lnf_1 + lnf_2 + (1 - lni_1)*(1 - lni_2)*(1-lnf_1)*(1-lnf_2);
    randnumber=scale*rand;
    
    if randnumber < lnf_1
        %wage_draw1(iter,time)=gendist(ff_1',1,1);
        urn = random('Beta', alphaf, betaf);
        transf_wage = wf_1(1) + urn*(wf_1(50)-wf_1(1));
        wage_aux=find(abs(wf_1 - transf_wage)<=min(abs(wf_1 - transf_wage)));
        wage_draw2(iter,time)=1;
        value(iter,time)=Wfn(wage_aux)*(Wfn(wage_aux)>Wnn)+Wnn*(Wfn(wage_aux)<=Wnn);
        state(iter,time)=3*(Wfn(wage_aux)>Wnn)+9*(Wfn(wage_aux)<=Wnn);
        wage_draw1(iter,time) = wage_aux*(Wfn(wage_aux)>Wnn) + (Wfn(wage_aux)<=Wnn);
        
    elseif randnumber >=lnf_1 && randnumber<(lni_1+lnf_1)
        %wage_draw1(iter,time)=gendist(fi_1',1,1);
        urn = random ('Beta', alphai, betai);
        transf_wage = wi_1(1) + urn*(wi_1(50)-wi_1(1));
        wage_aux=find(abs(wi_1 - transf_wage)<=min(abs(wi_1 - transf_wage)));
        wage_draw2(iter,time)=1;
        value(iter,time)=Win(wage_aux)*(Win(wage_aux)>Wnn)+Wnn*(Win(wage_aux)<=Wnn);
        state(iter,time)=7*(Win(wage_aux)>Wnn)+9*(Win(wage_aux)<=Wnn);
        wage_draw1(iter,time) = wage_aux*(Win(wage_aux)>Wnn) + (Win(wage_aux)<=Wnn);
        
    elseif randnumber >=(lni_1+lnf_1) && randnumber<(lnf_2+lni_1+lnf_1)
        wage_draw1(iter,time)=1;
        %wage_draw2(iter,time)=gendist(ff_2',1,1);
        urn = random ('Beta', alphaf, betaf);
        transf_wage = wf_2(1) + urn*(wf_2(50)-wf_2(1));
        wage_aux=find(abs(wf_2 - transf_wage)<=min(abs(wf_2 - transf_wage)));
        value(iter,time)=Wnf(wage_aux)*(Wnf(wage_aux)>Wnn)+Wnn*(Wnf(wage_aux)<=Wnn);
        state(iter,time)=5*(Wnf(wage_aux)>Wnn)+9*(Wnf(wage_aux)<=Wnn);
        wage_draw2(iter,time) = wage_aux*(Wnf(wage_aux)>Wnn) + (Wnf(wage_aux)<=Wnn);
        
    elseif randnumber >=(lnf_2+lni_1+lnf_1) && randnumber<(lni_2+lnf_2+lni_1+lnf_1)
        wage_draw1(iter,time)=1;
        urn = random ('Beta', alphai, betai);
        transf_wage = wi_2(1) + urn*(wi_2(50)-wi_2(1));
        wage_aux =find(abs(wi_2 - transf_wage)<=min(abs(wi_2 - transf_wage)));
        %wage_draw2(iter,time)=gendist(fi_2',1,1);
        value(iter,time)=Wni(wage_aux)*(Wni(wage_aux)>=Wnn)+Wnn*(Wni(wage_aux)<Wnn);
        state(iter,time)=8*(Wni(wage_aux)>=Wnn)+9*(Wni(wage_aux)<Wnn);
        wage_draw2(iter,time) = wage_aux*(Wni(wage_aux)>=Wnn)+(Wni(wage_aux)<Wnn);
    else
        wage_draw1(iter,time)=1;
        wage_draw2(iter,time)=1;
        value(iter,time)=Wnn;
        state(iter,time)=9;
        
    end
    
    
    time=2;
    
    while time<=maxtime
        
        
        
        
        %% FF in t-1
        
        
        
        if (state(iter,time-1) == 1)
            scale = df_1 + df_2 + lfi_1 + lfi_2 + (1 - df_1)*(1-df_2)*(1 - lfi_1)*(1 - lfi_2);
            randnumber=scale*rand;
            if  (randnumber < df_1)
                wage_draw1(iter,time)=wage_draw1(iter,time-1);
                wage_draw2(iter,time)=wage_draw2(iter,time-1);
                value(iter,time)=Wnf(wage_draw2(iter,time-1));
                state(iter,time)=5;
            end
            if ( randnumber >= df_1 && randnumber<(df_1+df_2))
                wage_draw1(iter,time)=wage_draw1(iter,time-1);
                wage_draw2(iter,time)=wage_draw2(iter,time-1);
                value(iter,time)=Wfn(wage_draw1(iter,time-1));
                state(iter,time)=3;
            end
            if ( randnumber >= (df_1+df_2) && randnumber<(lfi_1+df_1+df_2))
                %wage_draw1(iter,time)=gendist(fi_1',1,1);
                urn = random ('Beta', alphai, betai);
                transf_wage = wi_1(1) + urn*(wi_1(50)-wi_1(1));
                wage_aux=find(abs(wi_1 - transf_wage)<=min(abs(wi_1 - transf_wage)));
                wage_draw2(iter,time)=wage_draw2(iter,time-1);
                value(iter,time)=Wif(wage_aux,wage_draw2(iter,time-1));
                value(iter,time)=value(iter,time)*(value(iter,time)>value(iter,time-1))+value(iter,time-1)*(value(iter,time)<=value(iter,time-1));
                state(iter,time)=4*(value(iter,time)> value(iter,time-1))+state(iter,time-1)*(value(iter,time)<= value(iter,time-1));
                if state (iter,time) == state(iter,time-1)
                    wage_draw1(iter,time) = wage_draw1(iter,time - 1);
                else
                    wage_draw1(iter,time) = wage_aux;
                end
            end
            
            if (randnumber >= (lfi_1+df_1+df_2) && randnumber<(lfi_2+lfi_1+df_1+df_2))
                wage_draw1(iter,time)=wage_draw1(iter,time-1);
                %wage_draw2(iter,time)=gendist(fi_2',1,1);
                urn = random ('Beta', alphai, betai);
                transf_wage = wi_2(1) + urn*(wi_2(50)-wi_2(1));
                wage_aux=find(abs(wi_2 - transf_wage)<=min(abs(wi_2 - transf_wage)));
                value(iter,time)=Wfi(wage_draw1(iter,time-1),wage_aux);
                value(iter,time)=value(iter,time)*(value(iter,time)>value(iter,time-1))+value(iter,time-1)*(value(iter,time)<=value(iter,time-1));
                state(iter,time)=2*(value(iter,time)> value(iter,time-1))+state(iter,time-1)*(value(iter,time)<= value(iter,time-1));
                if state (iter,time) == state(iter,time-1)
                    wage_draw2(iter,time) = wage_draw2(iter,time - 1);
                else
                    wage_draw2(iter,time) = wage_aux;
                end
            end
            if (randnumber>=(lfi_2+lfi_1+df_1+df_2))
                wage_draw1(iter,time)=wage_draw1(iter,time-1);
                wage_draw2(iter,time)=wage_draw2(iter,time-1);
                value(iter,time)=value(iter,time-1);
                state(iter,time)=state(iter,time-1);
            end
        end
        %% FI em t-1
        if (state(iter,time-1)==2)
            scale = df_1 + di_2 + lfi_1 + lif_2 + (1 - di_2)*(1 -df_1)*(1-lfi_1)*(1 - lif_2);
            randnumber=scale*rand;
            
            if (randnumber < df_1)
                wage_draw1(iter,time)=wage_draw1(iter,time-1);
                wage_draw2(iter,time)=wage_draw2(iter,time-1);
                value(iter,time)=Wni(wage_draw2(iter,time-1));
                state(iter,time)=8;
            end
            if (randnumber >= df_1 && randnumber<(df_1+di_2))
                wage_draw1(iter,time)=wage_draw1(iter,time-1);
                wage_draw2(iter,time)=wage_draw2(iter,time-1);
                value(iter,time)=Wfn(wage_draw2(iter,time-1));
                state(iter,time)=3;
            end
            if randnumber >= (df_1+di_2) && randnumber<(lfi_1+df_1+di_2)
                %wage_draw1(iter,time)=gendist(fi_1',1,1);
                urn = random ('Beta', alphai, betai);
                transf_wage = wi_1(1) + urn*(wi_1(50)-wi_1(1));
                wage_aux=find(abs(wi_1 - transf_wage)<=min(abs(wi_1 - transf_wage)));
                wage_draw2(iter,time)=wage_draw2(iter,time-1);
                value(iter,time)=Wii(wage_aux,wage_draw2(iter,time-1));
                value(iter,time)=value(iter,time)*(value(iter,time)>value(iter,time-1))+value(iter,time-1)*(value(iter,time)<=value(iter,time-1));
                state(iter,time)=6*(value(iter,time)> value(iter,time-1))+state(iter,time-1)*(value(iter,time)<= value(iter,time-1));
                if state (iter,time) == state(iter,time-1)
                    wage_draw1(iter,time) = wage_draw1(iter,time - 1);
                else
                    wage_draw1(iter,time) = wage_aux;
                end
            end
            if randnumber >=(lfi_1+df_1+di_2) && randnumber<(lif_2+lfi_1+df_1+di_2)
                wage_draw1(iter,time)=wage_draw1(iter,time-1);
                %wage_draw2(iter,time)=gendist(ff_2',1,1);
                urn = random ('Beta', alphaf, betaf);
                transf_wage = wf_2(1) + urn*(wf_2(50)-wf_2(1));
                wage_aux=find(abs(wf_2 - transf_wage)<=min(abs(wf_2 - transf_wage)));
                value(iter,time)=Wff(wage_draw1(iter,time-1),wage_aux);
                value(iter,time)=value(iter,time)*(value(iter,time)>value(iter,time-1))+value(iter,time-1)*(value(iter,time)<=value(iter,time-1));
                state(iter,time)=1*(value(iter,time)> value(iter,time-1))+state(iter,time-1)*(value(iter,time)<= value(iter,time-1));
                if state (iter,time) == state(iter,time-1)
                    wage_draw2(iter,time) = wage_draw2(iter,time - 1);
                else
                    wage_draw2(iter,time) = wage_aux;
                end
            end
            if randnumber>=(lif_2+lfi_1+df_1+di_2)
                wage_draw1(iter,time)=wage_draw1(iter,time-1);
                wage_draw2(iter,time)=wage_draw2(iter,time-1);
                value(iter,time)=value(iter,time-1);
                state(iter,time)=state(iter,time-1);
                
            end
        end
        
        %% FN em t-1
        if state(iter,time-1)==3
            scale = df_1 + lfi_1 + lni_2 +lnf_2 + (1-df_1)*(1-lfi_1)*(1-lni_2)*(1-lnf_2);
            randnumber=scale*rand;
            
            if randnumber < df_1*(1-p_2)
                wage_draw1(iter,time)=wage_draw1(iter,time-1);
                wage_draw2(iter,time)=wage_draw2(iter,time-1);
                value(iter,time)=Wnn;
                state(iter,time)=9;
            end
            if randnumber >= df_1*(1-p_2) && randnumber < df_1
                wage_draw1(iter,time)=wage_draw1(iter,time-1);
                %wage_draw2(iter,time)=gendist(fi_2',1,1);
                urn = random ('Beta', alphai, betai);
                transf_wage = wi_2(1) + urn*(wi_2(50)-wi_2(1));
                wage_aux=find(abs(wi_2 - transf_wage)<=min(abs(wi_2 - transf_wage)));
                value(iter,time)=Wni(wage_aux);
                value(iter,time)=value(iter,time)*(value(iter,time)>Wnn)+Wnn*(value(iter,time)<=Wnn);
                state(iter,time)=8*(value(iter,time)> Wnn)+9*(value(iter,time)<=Wnn);
                if state (iter,time) == 9
                    wage_draw1(iter,time) = wage_draw1(iter,time - 1);
                else
                    wage_draw1(iter,time) = wage_aux;
                end
            end
            if randnumber >= df_1 && randnumber<(lfi_1+df_1)
                %wage_draw1(iter,time)=gendist(fi_1',1,1);
                urn = random ('Beta', alphai, betai);
                transf_wage = wi_1(1) + urn*(wi_1(50)-wi_1(1));
                wage_aux=find(abs(wi_1 - transf_wage)<=min(abs(wi_1 - transf_wage)));
                wage_draw2(iter,time)=wage_draw2(iter,time-1);
                value(iter,time)=Win(wage_aux);
                value(iter,time)=value(iter,time)*(value(iter,time)>value(iter,time-1))+value(iter,time-1)*(value(iter,time)<=value(iter,time-1));
                state(iter,time)=7*(value(iter,time)> value(iter,time-1))+state(iter,time-1)*(value(iter,time)<= value(iter,time-1));  
                if state (iter,time) == state(iter,time-1)
                    wage_draw1(iter,time) = wage_draw1(iter,time - 1);
                else
                    wage_draw1(iter,time) = wage_aux;
                end
            end
            if randnumber >= (lfi_1+df_1) && randnumber<(lnf_2+lfi_1+df_1)
                wage_draw1(iter,time)=wage_draw1(iter,time-1);
                % wage_draw2(iter,time)=gendist(ff_2',1,1);
                urn = random ('Beta', alphaf, betaf);
                transf_wage = wf_2(1) + urn*(wf_2(50)-wf_2(1));
                wage_aux=find(abs(wf_2 - transf_wage)<=min(abs(wf_2 - transf_wage)));
                value(iter,time)=Wff(wage_draw1(iter,time-1),wage_aux);
                value(iter,time)=value(iter,time)*(value(iter,time)>value(iter,time-1))+value(iter,time-1)*(value(iter,time)<=value(iter,time-1));
                state(iter,time)=1*(value(iter,time)> value(iter,time-1))+state(iter,time-1)*(value(iter,time)<= value(iter,time-1));
                if state (iter,time) == state(iter,time-1)
                    wage_draw2(iter,time) = wage_draw2(iter,time - 1);
                else
                    wage_draw2(iter,time) = wage_aux;
                end
            end
            if randnumber >= (lfi_1+df_1 + lnf_2) && randnumber<(lni_2+lnf_2+lfi_1+df_1)
                wage_draw1(iter,time)=wage_draw1(iter,time-1);
                urn = random ('Beta', alphai, betai);
                transf_wage = wi_2(1) + urn*(wi_2(50)-wi_2(1));
                wage_aux=find(abs(wi_2 - transf_wage)<=min(abs(wi_2 - transf_wage)));
                %wage_draw2(iter,time)=gendist(fi_2',1,1);
                value(iter,time)=Wfi(wage_draw1(iter,time-1),wage_aux);
                value(iter,time)=value(iter,time)*(value(iter,time)>value(iter,time-1))+value(iter,time-1)*(value(iter,time)<=value(iter,time-1));
                state(iter,time)=2*(value(iter,time)> value(iter,time-1))+state(iter,time-1)*(value(iter,time)<= value(iter,time-1));
                if state (iter,time) == state(iter,time-1)
                    wage_draw2(iter,time) = wage_draw2(iter,time - 1);
                else
                    wage_draw2(iter,time) = wage_aux;
                end
            end
            if randnumber >= (lni_2+lnf_2+lfi_1+df_1)
                wage_draw1(iter,time)=wage_draw1(iter,time-1);
                wage_draw2(iter,time)=wage_draw2(iter,time-1);
                value(iter,time)=value(iter,time-1);
                state(iter,time)=state(iter,time-1);
            end
        end
        %% IF em t-1
        if (state(iter,time-1)==4)
            scale = df_2 + di_1 + lfi_2 + lif_1 + (1 - di_1)*(1 -df_2)*(1-lfi_2)*(1 - lif_1);
            randnumber=scale*rand;
            
            if randnumber < di_1
                wage_draw1(iter,time)=wage_draw1(iter,time-1);
                wage_draw2(iter,time)=wage_draw2(iter,time-1);
                value(iter,time)=Wnf(wage_draw2(iter,time-1));
                state(iter,time)=5;
            end
            if randnumber >= di_1 && randnumber<(df_2+di_1)
                wage_draw1(iter,time)=wage_draw1(iter,time-1);
                wage_draw2(iter,time)=wage_draw2(iter,time-1);
                value(iter,time)=Win(wage_draw2(iter,time-1));
                state(iter,time)=7;
            end
            if randnumber >= (df_2+di_1) && randnumber<(lif_1+df_2+di_1)
                % wage_draw1(iter,time)=gendist(ff_1',1,1);
                urn = random('Beta', alphaf, betaf);
                transf_wage = wf_1(1) + urn*(wf_1(50)-wf_1(1));
                wage_aux=find(abs(wf_1 - transf_wage)<=min(abs(wf_1 - transf_wage)));
                wage_draw2(iter,time)=wage_draw2(iter,time-1);
                value(iter,time)=Wff(wage_aux,wage_draw2(iter,time-1));
                value(iter,time)=value(iter,time)*(value(iter,time)>value(iter,time-1))+value(iter,time-1)*(value(iter,time)<=value(iter,time-1));
                state(iter,time)=1*(value(iter,time)> value(iter,time-1))+state(iter,time-1)*(value(iter,time)<= value(iter,time-1));
                if state (iter,time) == state(iter,time-1)
                    wage_draw1(iter,time) = wage_draw1(iter,time - 1);
                else
                    wage_draw1(iter,time) = wage_aux;
                end
            end
            if randnumber >=(lif_1+df_2+di_1) && randnumber<(lfi_2+lif_1+df_2+di_1)
                wage_draw1(iter,time)=wage_draw1(iter,time-1);
                urn = random ('Beta', alphai, betai);
                transf_wage = wi_2(1) + urn*(wi_2(50)-wi_2(1));
                wage_aux=find(abs(wi_2 - transf_wage)<=min(abs(wi_2 - transf_wage)));
                %wage_draw2(iter,time)=gendist(fi_2',1,1);
                value(iter,time)=Wii(wage_draw1(iter,time-1),wage_aux);
                value(iter,time)=value(iter,time)*(value(iter,time)>value(iter,time-1))+value(iter,time-1)*(value(iter,time)<=value(iter,time-1));
                state(iter,time)=6*(value(iter,time)> value(iter,time-1))+state(iter,time-1)*(value(iter,time)<= value(iter,time-1));
                if state (iter,time) == state(iter,time-1)
                    wage_draw2(iter,time) = wage_draw2(iter,time - 1);
                else
                    wage_draw2(iter,time) = wage_aux;
                end
            end
            if randnumber>=(lfi_2+lif_1+df_2+di_1)
                wage_draw1(iter,time)=wage_draw1(iter,time-1);
                wage_draw2(iter,time)=wage_draw2(iter,time-1);
                value(iter,time)=value(iter,time-1);
                state(iter,time)=state(iter,time-1);
            end
        end
        %% NF em t-1
        if state(iter,time-1)==5
            scale = df_2 + lfi_2 + lnf_1 + lni_1 + (1 - lni_1)*(1 - lnf_1)*(1-lfi_2)*(1 - df_2);
            randnumber=scale*rand;
            
            if randnumber < df_2*(1-p_1)
                wage_draw1(iter,time)=wage_draw1(iter,time-1);
                wage_draw2(iter,time)=wage_draw2(iter,time-1);
                value(iter,time)=Wnn;
                state(iter,time)=9;
            end
            if randnumber >= df_2*(1-p_1) && randnumber < df_2
                %wage_draw1(iter,time)=gendist(fi_1',1,1);
                urn = random ('Beta', alphai, betai);
                transf_wage = wi_1(1) + urn*(wi_1(50)-wi_1(1));
                wage_aux=find(abs(wi_1 - transf_wage)<=min(abs(wi_1 - transf_wage)));
                wage_draw2(iter,time)=wage_draw2(iter,time-1);
                value(iter,time)=Win(wage_aux);
                value(iter,time)=value(iter,time)*(value(iter,time)>Wnn)+Wnn*(value(iter,time)<=Wnn);
                state(iter,time)=7*(value(iter,time)> Wnn)+9*(value(iter,time)<=Wnn);
                if state (iter,time) == 9
                    wage_draw1(iter,time) = wage_draw1(iter,time - 1);
                else
                    wage_draw1(iter,time) = wage_aux;
                end
            end
            if randnumber >= df_2 && randnumber<(lfi_2+df_2)
                wage_draw1(iter,time)=wage_draw1(iter,time-1);
                urn = random ('Beta', alphai, betai);
                transf_wage = wi_2(1) + urn*(wi_2(50)-wi_2(1));
                wage_aux=find(abs(wi_2 - transf_wage)<=min(abs(wi_2 - transf_wage)));
                %  wage_draw2(iter,time)=gendist(fi_2',1,1);
                value(iter,time)=Wni(wage_aux);
                value(iter,time)=value(iter,time)*(value(iter,time)>value(iter,time-1))+value(iter,time-1)*(value(iter,time)<=value(iter,time-1));
                state(iter,time)=8*(value(iter,time)> value(iter,time-1))+state(iter,time-1)*(value(iter,time)<= value(iter,time-1));
                   if state (iter,time) == state(iter,time-1)
                    wage_draw2(iter,time) = wage_draw2(iter,time - 1);
                else
                    wage_draw2(iter,time) = wage_aux;
                end
            end
            if randnumber >= (lfi_2+df_2) && randnumber<(lnf_1+lfi_2+df_2)
                urn = random('Beta', alphaf, betaf);
                transf_wage = wf_1(1) + urn*(wf_1(50)-wf_1(1));
                wage_aux=find(abs(wf_1 - transf_wage)<=min(abs(wf_1 - transf_wage)));
                %wage_draw1(iter,time)=gendist(ff_1',1,1);
                wage_draw2(iter,time)=wage_draw2(iter,time-1);
                value(iter,time)=Wff(wage_aux,wage_draw2(iter,time-1));
                value(iter,time)=value(iter,time)*(value(iter,time)>value(iter,time-1))+value(iter,time-1)*(value(iter,time)<=value(iter,time-1));
                state(iter,time)=1*(value(iter,time)> value(iter,time-1))+state(iter,time-1)*(value(iter,time)<= value(iter,time-1));
                if state (iter,time) == state(iter,time-1)
                    wage_draw1(iter,time) = wage_draw1(iter,time - 1);
                else
                    wage_draw1(iter,time) = wage_aux;
                end
            end
            if randnumber >= (lnf_1+lfi_2+df_2) && randnumber<(lni_1+lnf_1+lfi_2+df_2)
                % wage_draw1(iter,time)=gendist(fi_1',1,1);
                urn = random ('Beta', alphai, betai);
                transf_wage = wi_1(1) + urn*(wi_1(50)-wi_1(1));
                wage_aux=find(abs(wi_1 - transf_wage)<=min(abs(wi_1 - transf_wage)));
                wage_draw2(iter,time)=wage_draw2(iter,time-1);
                value(iter,time)=Wif(wage_aux,wage_draw2(iter,time));
                value(iter,time)=value(iter,time)*(value(iter,time)>value(iter,time-1))+value(iter,time-1)*(value(iter,time)<=value(iter,time-1));
                state(iter,time)=4*(value(iter,time)> value(iter,time-1))+state(iter,time-1)*(value(iter,time)<= value(iter,time-1));
                if state (iter,time) == state(iter,time-1)
                    wage_draw1(iter,time) = wage_draw1(iter,time - 1);
                else
                    wage_draw1(iter,time) = wage_aux;
                end
            end
            if randnumber >= (lni_1+lnf_1+lfi_2+df_2)
                wage_draw1(iter,time)=wage_draw1(iter,time-1);
                wage_draw2(iter,time)=wage_draw2(iter,time-1);
                value(iter,time)=value(iter,time-1);
                state(iter,time)=state(iter,time-1);
            end
        end
        %% II em t-1
        if state(iter,time-1)==6
            scale = di_1 + di_2 + lif_1 + lif_2 + (1-di_1)*(1-di_2)*(1-lif_1)*(1-lif_2);
            randnumber=scale*rand;
            
            if randnumber < di_1
                wage_draw1(iter,time)=wage_draw1(iter,time-1);
                wage_draw2(iter,time)=wage_draw2(iter,time-1);
                value(iter,time)=Wni(wage_draw2(iter,time-1));
                state(iter,time)=8;
            end
            if randnumber >= di_1 && randnumber<(di_1+di_2)
                wage_draw1(iter,time)=wage_draw1(iter,time-1);
                wage_draw2(iter,time)=wage_draw2(iter,time-1);
                value(iter,time)=Win(wage_draw1(iter,time-1));
                state(iter,time)=7;
            end
            if randnumber >= (di_1+di_2) && randnumber<(lif_1+di_1+di_2)
                urn = random('Beta', alphaf, betaf);
                transf_wage = wf_1(1) + urn*(wf_1(50)-wf_1(1));
                wage_aux=find(abs(wf_1 - transf_wage)<=min(abs(wf_1 - transf_wage)));
                % wage_draw1(iter,time)=gendist(ff_1',1,1);
                wage_draw2(iter,time)=wage_draw2(iter,time-1);
                value(iter,time)=Wfi(wage_aux,wage_draw2(iter,time-1));
                value(iter,time)=value(iter,time)*(value(iter,time)>value(iter,time-1))+value(iter,time-1)*(value(iter,time)<=value(iter,time-1));
                state(iter,time)=2*(value(iter,time)> value(iter,time-1))+state(iter,time-1)*(value(iter,time)<= value(iter,time-1));
                if state (iter,time) == state(iter,time-1)
                    wage_draw1(iter,time) = wage_draw1(iter,time - 1);
                else
                    wage_draw1(iter,time) = wage_aux;
                end
            end
            if randnumber >= (lif_1+di_1+di_2) && randnumber<(lif_2+lif_1+di_1+di_2)
                wage_draw1(iter,time)=wage_draw1(iter,time-1);
                urn = random('Beta', alphaf, betaf);
                transf_wage = wf_2(1) + urn*(wf_2(50)-wf_2(1));
                wage_aux=find(abs(wf_2 - transf_wage)<=min(abs(wf_2 - transf_wage)));
                value(iter,time)=Wif(wage_draw1(iter,time-1),wage_aux);
                value(iter,time)=value(iter,time)*(value(iter,time)>value(iter,time-1))+value(iter,time-1)*(value(iter,time)<=value(iter,time-1));
                state(iter,time)=4*(value(iter,time)> value(iter,time-1))+state(iter,time-1)*(value(iter,time)<= value(iter,time-1));
                if state (iter,time) == state(iter,time-1)
                    wage_draw2(iter,time) = wage_draw2(iter,time - 1);
                else
                    wage_draw2(iter,time) = wage_aux;
                end
            end
            if randnumber>=(lif_2+lif_1+di_1+di_2)
                wage_draw1(iter,time)=wage_draw1(iter,time-1);
                wage_draw2(iter,time)=wage_draw2(iter,time-1);
                value(iter,time)=value(iter,time-1);
                state(iter,time)=state(iter,time-1);
            end
        end
        %% IN em t-1
        if state(iter,time-1)==7
            scale = di_1 + lif_1 + lnf_2 + lni_2 + (1 - di_1)*(1 - lif_1)*(1 - lnf_2)*(1- lni_1);
            randnumber=scale*rand;
            
            if randnumber < di_1*(1-q_2)
                wage_draw1(iter,time)=wage_draw1(iter,time-1);
                wage_draw2(iter,time)=wage_draw2(iter,time-1);
                value(iter,time)=Wnn;
                state(iter,time)=9;
            end
            if randnumber >= di_1*(1-q_2) && randnumber < di_1
                wage_draw1(iter,time)=wage_draw1(iter,time-1);
                urn = random ('Beta', alphai, betai);
                transf_wage = wi_2(1) + urn*(wi_2(50)-wi_2(1));
                wage_aux=find(abs(wi_2 - transf_wage)<=min(abs(wi_2 - transf_wage)));
                % wage_draw2(iter,time)=gendist(fi_2',1,1);
                value(iter,time)=Wni(wage_aux);
                value(iter,time)=value(iter,time)*(value(iter,time)>Wnn)+Wnn*(value(iter,time)<=Wnn);
                state(iter,time)=8*(value(iter,time)> Wnn)+9*(value(iter,time)<=Wnn);
                if state (iter,time) == 9
                    wage_draw2(iter,time) = wage_draw2(iter,time - 1);
                else
                    wage_draw2(iter,time) = wage_aux;
                end
            end
            if randnumber >= di_1 && randnumber<(lif_1+di_1)
                urn = random('Beta', alphaf, betaf);
                transf_wage = wf_1(1) + urn*(wf_1(50)-wf_1(1));
                wage_aux=find(abs(wf_1 - transf_wage)<=min(abs(wf_1 - transf_wage)));
                %wage_draw1(iter,time)=gendist(ff_1',1,1);
                wage_draw2(iter,time)=wage_draw2(iter,time-1);
                value(iter,time)=Wfn(wage_aux);
                value(iter,time)=value(iter,time)*(value(iter,time)>value(iter,time-1))+value(iter,time-1)*(value(iter,time)<=value(iter,time-1));
                state(iter,time)=3*(value(iter,time)> value(iter,time-1))+state(iter,time-1)*(value(iter,time)<= value(iter,time-1));
                if state (iter,time) == state(iter,time-1)
                    wage_draw1(iter,time) = wage_draw1(iter,time - 1);
                else
                    wage_draw1(iter,time) = wage_aux;
                end
            end
            if randnumber >= (lif_1+di_1) && randnumber<(lnf_2+lif_1+di_1)
                wage_draw1(iter,time)=wage_draw1(iter,time-1);
                urn = random ('Beta', alphaf, betaf);
                transf_wage = wf_2(1) + urn*(wf_2(50)-wf_2(1));
                wage_aux=find(abs(wf_2 - transf_wage)<=min(abs(wf_2 - transf_wage)));
                %wage_draw2(iter,time)=gendist(ff_2',1,1);
                value(iter,time)=Wif(wage_draw1(iter,time-1),wage_aux);
                value(iter,time)=value(iter,time)*(value(iter,time)>value(iter,time-1))+value(iter,time-1)*(value(iter,time)<=value(iter,time-1));
                state(iter,time)=4*(value(iter,time)> value(iter,time-1))+state(iter,time-1)*(value(iter,time)<= value(iter,time-1));
                if state (iter,time) == state(iter,time-1)
                    wage_draw2(iter,time) = wage_draw2(iter,time - 1);
                else
                    wage_draw2(iter,time) = wage_aux;
                end
            end
            if randnumber >= (lnf_2+lif_1+di_1) && randnumber<(lni_2+lnf_2+lif_1+di_1)
                wage_draw1(iter,time)=wage_draw1(iter,time-1);
                urn = random ('Beta', alphai, betai);
                transf_wage = wi_2(1) + urn*(wi_2(50)-wi_2(1));
                wage_aux=find(abs(wi_2 - transf_wage)<=min(abs(wi_2 - transf_wage)));
                %wage_draw2(iter,time)=gendist(fi_2',1,1);
                value(iter,time)=Wii(wage_draw1(iter,time-1),wage_aux);
                value(iter,time)=value(iter,time)*(value(iter,time)>value(iter,time-1))+value(iter,time-1)*(value(iter,time)<=value(iter,time-1));
                state(iter,time)=6*(value(iter,time)> value(iter,time-1))+state(iter,time-1)*(value(iter,time)<= value(iter,time-1));
                if state (iter,time) == state(iter,time-1)
                    wage_draw2(iter,time) = wage_draw2(iter,time - 1);
                else
                    wage_draw2(iter,time) = wage_aux;
                end
            end
            if randnumber >= (lni_2+lnf_2+lif_1+di_1)
                wage_draw1(iter,time)=wage_draw1(iter,time-1);
                wage_draw2(iter,time)=wage_draw2(iter,time-1);
                value(iter,time)=value(iter,time-1);
                state(iter,time)=state(iter,time-1);
            end
        end
        %% NI em t-1
        if state(iter,time-1)==8
            scale = di_2 + lif_2 + lni_1 + lnf_1 + (1 - di_2)*(1 - lif_2)*(1 - lni_1)*(1 - lnf_1);
            randnumber=scale*rand;
            
            if randnumber < di_2*(1-q_1)
                wage_draw1(iter,time)=wage_draw1(iter,time-1);
                wage_draw2(iter,time)=wage_draw2(iter,time-1);
                value(iter,time)=Wnn;
                state(iter,time)=9;
            end
            if randnumber >= di_2*(1-q_1) && randnumber < di_2
                %wage_draw1(iter,time)=gendist(fi_1',1,1);
                urn = random ('Beta', alphai, betai);
                transf_wage = wi_1(1) + urn*(wi_1(50)-wi_1(1));
                wage_aux=find(abs(wi_1 - transf_wage)<=min(abs(wi_1 - transf_wage)));
                wage_draw2(iter,time)=wage_draw2(iter,time-1);
                value(iter,time)=Win(wage_aux);
                value(iter,time)=value(iter,time)*(value(iter,time)>Wnn)+Wnn*(value(iter,time)<=Wnn);
                state(iter,time)=7*(value(iter,time)> Wnn)+9*(value(iter,time)<=Wnn);
                if state (iter,time) == 9
                    wage_draw1(iter,time) = wage_draw1(iter,time - 1);
                else
                    wage_draw1(iter,time) = wage_aux;
                end
            end
            if randnumber >= di_2 && randnumber<(lif_2+di_2)
                wage_draw1(iter,time)=wage_draw1(iter,time-1);
                %wage_draw2(iter,time)=gendist(ff_2',1,1);
                urn = random ('Beta', alphaf, betaf);
                transf_wage = wf_2(1) + urn*(wf_2(50)-wf_2(1));
                wage_aux=find(abs(wf_2 - transf_wage)<=min(abs(wf_2 - transf_wage)));
                value(iter,time)=Wnf(wage_aux);
                value(iter,time)=value(iter,time)*(value(iter,time)>value(iter,time-1))+value(iter,time-1)*(value(iter,time)<=value(iter,time-1));
                state(iter,time)=5*(value(iter,time)> value(iter,time-1))+state(iter,time-1)*(value(iter,time)<= value(iter,time-1));
                 if state (iter,time) == state(iter,time-1)
                    wage_draw2(iter,time) = wage_draw2(iter,time - 1);
                else
                    wage_draw2(iter,time) = wage_aux;
                end
            end
            if randnumber >= (lif_2+di_2) && randnumber<(lnf_1+lif_2+di_2)
                urn = random('Beta', alphaf, betaf);
                transf_wage = wf_1(1) + urn*(wf_1(50)-wf_1(1));
                wage_aux=find(abs(wf_1 - transf_wage)<=min(abs(wf_1 - transf_wage)));
                % wage_draw1(iter,time)=gendist(ff_1',1,1);
                wage_draw2(iter,time)=wage_draw2(iter,time-1);
                value(iter,time)=Wfi(wage_aux,wage_draw2(iter,time));
                value(iter,time)=value(iter,time)*(value(iter,time)>value(iter,time-1))+value(iter,time-1)*(value(iter,time)<=value(iter,time-1));
                state(iter,time)=2*(value(iter,time)> value(iter,time-1))+state(iter,time-1)*(value(iter,time)<= value(iter,time-1));
                 if state (iter,time) == state(iter,time-1)
                    wage_draw1(iter,time) = wage_draw1(iter,time - 1);
                else
                    wage_draw1(iter,time) = wage_aux;
                end
            end
            if randnumber >= (lnf_1+lif_2+di_2) && randnumber<(lni_1+lnf_1+lif_2+di_2)
                %wage_draw1(iter,time)=gendist(fi_1',1,1);
                urn = random ('Beta', alphai, betai);
                transf_wage = wi_1(1) + urn*(wi_1(50)-wi_1(1));
                wage_aux=find(abs(wi_1 - transf_wage)<=min(abs(wi_1 - transf_wage)));
                wage_draw2(iter,time)=wage_draw2(iter,time-1);
                value(iter,time)=Wii(wage_aux,wage_draw2(iter,time));
                value(iter,time)=value(iter,time)*(value(iter,time)>value(iter,time-1))+value(iter,time-1)*(value(iter,time)<=value(iter,time-1));
                state(iter,time)=6*(value(iter,time)> value(iter,time-1))+state(iter,time-1)*(value(iter,time)<= value(iter,time-1));
                 if state (iter,time) == state(iter,time-1)
                    wage_draw1(iter,time) = wage_draw1(iter,time - 1);
                else
                    wage_draw1(iter,time) = wage_aux;
                end
            end
            if randnumber >= (lni_1+lnf_1+lif_2+di_2)
                wage_draw1(iter,time)=wage_draw1(iter,time-1);
                wage_draw2(iter,time)=wage_draw2(iter,time-1);
                value(iter,time)=value(iter,time-1);
                state(iter,time)=state(iter,time-1);
            end
        end
        % NN in t-1
        if state(iter,time-1) ==9
            scale = lnf_1 + lni_1 + lni_2 + lnf_2 + (1 - lni_2)*(1 - lni_1)*(1 - lnf_2)*(1 - lnf_1);
            randnumber=scale*rand;
            
            if randnumber < lnf_1
                %wage_draw1(iter,time)=gendist(ff_1',1,1);
                urn = random('Beta', alphaf, betaf);
                transf_wage = wf_1(1) + urn*(wf_1(50)-wf_1(1));
                wage_aux=find(abs(wf_1 - transf_wage)<=min(abs(wf_1 - transf_wage)));
                wage_draw2(iter,time)=1;
                value(iter,time)=Wfn(wage_aux)*(Wfn(wage_aux)>Wnn)+Wnn*(Wfn(wage_aux)<=Wnn);
                state(iter,time)=3*(Wfn(wage_aux)>Wnn)+9*(Wfn(wage_aux)<=Wnn);
                 if state (iter,time) == state(iter,time-1)
                    wage_draw1(iter,time) = wage_draw1(iter,time - 1);
                else
                    wage_draw1(iter,time) = wage_aux;
                end
            end
            if randnumber >=lnf_1 && randnumber<(lni_1+lnf_1)
                %wage_draw1(iter,time)=gendist(fi_1',1,1);
                urn = random ('Beta', alphai, betai);
                transf_wage = wi_1(1) + urn*(wi_1(50)-wi_1(1));
                wage_aux=find(abs(wi_1 - transf_wage)<=min(abs(wi_1 - transf_wage)));
                wage_draw2(iter,time)=1;
                value(iter,time)=Win(wage_aux)*(Win(wage_aux)>Wnn)+Wnn*(Win(wage_aux)<=Wnn);
                state(iter,time)=7*(Win(wage_aux)>Wnn)+9*(Win(wage_aux)<=Wnn);
                 if state (iter,time) == state(iter,time-1)
                    wage_draw1(iter,time) = wage_draw1(iter,time - 1);
                else
                    wage_draw1(iter,time) = wage_aux;
                end
            end
            if randnumber >=(lni_1+lnf_1) && randnumber<(lnf_2+lni_1+lnf_1)
                wage_draw1(iter,time)=1;
                %wage_draw2(iter,time)=gendist(ff_2',1,1);
                urn = random ('Beta', alphaf, betaf);
                transf_wage = wf_2(1) + urn*(wf_2(50)-wf_2(1));
                wage_aux=find(abs(wf_2 - transf_wage)<=min(abs(wf_2 - transf_wage)));
                value(iter,time)=Wnf(wage_aux)*(Wnf(wage_aux)>Wnn)+Wnn*(Wnf(wage_aux)<=Wnn);
                state(iter,time)=5*(Wnf(wage_aux)>Wnn)+9*(Wnf(wage_aux)<=Wnn);
                 if state (iter,time) == state(iter,time-1)
                    wage_draw2(iter,time) = wage_draw2(iter,time - 1);
                else
                    wage_draw2(iter,time) = wage_aux;
                end
            end
            if randnumber >=(lnf_2+lni_1+lnf_1) && randnumber<(lni_2+lnf_2+lni_1+lnf_1)
                wage_draw1(iter,time)=1;
                urn = random ('Beta', alphai, betai);
                transf_wage = wi_2(1) + urn*(wi_2(50)-wi_2(1));
                wage_aux=find(abs(wi_2 - transf_wage)<=min(abs(wi_2 - transf_wage)));
                %  wage_draw2(iter,time)=gendist(fi_2',1,1);
                value(iter,time)=Wni(wage_aux)*(Wni(wage_aux)>=Wnn)+Wnn*(Wni(wage_aux)<Wnn);
                state(iter,time)=8*(Wni(wage_aux)>=Wnn)+9*(Wni(wage_aux)<Wnn);
                  if state (iter,time) == state(iter,time-1)
                    wage_draw2(iter,time) = wage_draw2(iter,time - 1);
                else
                    wage_draw2(iter,time) = wage_aux;
                end
            end
            if randnumber>=(lni_2+lnf_2+lni_1+lnf_1)
                value(iter,time)=Wnn;
                state(iter,time)=9;
            end
        end
        
        
        time=time+1;
        
    end
    
end

% Moments calculation

% use only last 2 quarters of simulation for each household
FF=sum(state(:,end-1)==1)/N_households;
FI=sum(state(:,end-1)==2)/N_households;
FN=sum(state(:,end-1)==3)/N_households;
IF=sum(state(:,end-1)==4)/N_households;
NF=sum(state(:,end-1)==5)/N_households;
II=sum(state(:,end-1)==6)/N_households;
IN=sum(state(:,end-1)==7)/N_households;
NI=sum(state(:,end-1)==8)/N_households;
NN=sum(state(:,end-1)==9)/N_households;

HHstates=[FF;FI;FN;IF;NF;II;IN;NI;NN]; % joint-job state

hh_job_status1=(state(:,end-1)>=6); % save data to estimate gamma
save hh_job_status1 hh_job_status1

wage_formal_1=wf_1(wage_draw1(:,end-1));
wage_informal_1=wi_1(wage_draw1(:,end-1));
wage_formal_2=wf_2(wage_draw2(:,end-1));
wage_informal_2=wi_2(wage_draw2(:,end-1));

indicewf_1=find(state(:,end-1)<=3);
wage_formal_1=sort(wage_formal_1(indicewf_1));
indicewi_1=find(state(:,end-1)==4|state(:,end-1)==6|state(:,end-1)==7);
wage_informal_1=sort(wage_informal_1(indicewi_1));
indicewf_2=find(state(:,end-1)==1|state(:,end-1)==4|state(:,end-1)==5);
wage_formal_2=sort(wage_formal_2(indicewf_2));
indicewi_2=find(state(:,end-1)==2|state(:,end-1)==6|state(:,end-1)==8);
wage_informal_2=sort(wage_informal_2(indicewi_2));

sdwi_1 = std(wage_informal_1);
sdwi_2 = std(wage_informal_2);
sdwf_1 = std(wage_formal_1);
sdwf_2 = std(wage_formal_2);

meanwf_1=mean(wage_formal_1);

p10wf_1=prctile(wage_formal_1,10);
p25wf_1=prctile(wage_formal_1,25);
p50wf_1=prctile(wage_formal_1,50);
p75wf_1=prctile(wage_formal_1,75);
p90wf_1=prctile(wage_formal_1,90);



meanwf_2=mean(wage_formal_2);

p10wf_2=prctile(wage_formal_2,10);
p25wf_2=prctile(wage_formal_2,25);
p50wf_2=prctile(wage_formal_2,50);
p75wf_2=prctile(wage_formal_2,75);
p90wf_2=prctile(wage_formal_2,90);

meanwi_1=mean(wage_informal_1);

p10wi_1=prctile(wage_informal_1,10);
p25wi_1=prctile(wage_informal_1,25);
p50wi_1=prctile(wage_informal_1,50);
p75wi_1=prctile(wage_informal_1,75);
p90wi_1=prctile(wage_informal_1,90);


meanwi_2=mean(wage_informal_2);

p10wi_2=prctile(wage_informal_2,10);
p25wi_2=prctile(wage_informal_2,25);
p50wi_2=prctile(wage_informal_2,50);
p75wi_2=prctile(wage_informal_2,75);
p90wi_2=prctile(wage_informal_2,90);


wage_moments=[meanwi_1; meanwf_1; meanwi_2 ; meanwf_2; sdwi_1; sdwf_1; sdwi_2; sdwf_2; p10wi_1 ; p10wf_1; p10wi_2 ; p10wf_2;...
    p25wi_1; p25wf_1; p25wi_2 ; p25wf_2; p50wi_1; p50wf_1; p50wi_2 ; p50wf_2; p75wi_1; p75wf_1; p75wi_2 ; p75wf_2;...
    p90wi_1; p90wf_1 ; p90wi_2 ; p90wf_2]; % wages

% conditional transitions
formal_1_time1=(state(:,end-1)==1|state(:,end-1)==2|state(:,end-1)==3);
informal_1_time1=(state(:,end-1)==4|state(:,end-1)==6|state(:,end-1)==7);
nonemp_1_time1=(state(:,end-1)==5|state(:,end-1)==8|state(:,end-1)==9);
formal_2_time1=(state(:,end-1)==1|state(:,end-1)==4|state(:,end-1)==5);
informal_2_time1=(state(:,end-1)==2|state(:,end-1)==6|state(:,end-1)==8);
nonemp_2_time1=(state(:,end-1)==3|state(:,end-1)==7|state(:,end-1)==9);

formal_1_time2=(state(:,end)==1|state(:,end)==2|state(:,end)==3);
informal_1_time2=(state(:,end)==4|state(:,end)==6|state(:,end)==7);
nonemp_1_time2=(state(:,end)==5|state(:,end)==8|state(:,end)==9);
formal_2_time2=(state(:,end)==1|state(:,end)==4|state(:,end)==5);
informal_2_time2=(state(:,end)==2|state(:,end)==6|state(:,end)==8);
nonemp_2_time2=(state(:,end)==3|state(:,end)==7|state(:,end)==9);

Dfn_1=sum(formal_1_time1==1 & nonemp_1_time2==1)/sum(formal_1_time1==1);
Din_1=sum(informal_1_time1==1 & nonemp_1_time2==1)/sum(informal_1_time1==1);
Dnf_1=sum(nonemp_1_time1==1 & formal_1_time2==1)/sum(nonemp_1_time1==1);
Dni_1=sum(nonemp_1_time1==1 & informal_1_time2==1)/sum(nonemp_1_time1==1);
Dfi_1=sum(formal_1_time1==1 & informal_1_time2==1)/sum(formal_1_time1==1);
Dif_1=sum(informal_1_time1==1 & formal_1_time2==1)/sum(informal_1_time1==1);
Dfn_2=sum(formal_2_time1==1 & nonemp_2_time2==1)/sum(formal_2_time1==1);
Din_2=sum(informal_2_time1==1 & nonemp_2_time2==1)/sum(informal_2_time1==1);
Dnf_2=sum(nonemp_2_time1==1 & formal_2_time2==1)/sum(nonemp_2_time1==1);
Dni_2=sum(nonemp_2_time1==1 & informal_2_time2==1)/sum(nonemp_2_time1==1);
Dfi_2=sum(formal_2_time1==1 & informal_2_time2==1)/sum(formal_2_time1==1);
Dif_2=sum(informal_2_time1==1 & formal_2_time2==1)/sum(informal_2_time1==1);
Dni_1_Dfn_2=sum(nonemp_1_time1==1 & informal_1_time2==1 & formal_2_time1==1 & nonemp_2_time2==1)/sum(nonemp_1_time1==1);
Dni_1_Din_2=sum(nonemp_1_time1==1 & informal_1_time2==1 & informal_2_time1==1 & nonemp_2_time2==1)/sum(nonemp_1_time1==1);
Dni_2_Dfn_1=sum(nonemp_2_time1==1 & informal_2_time2==1 & formal_1_time1==1 & nonemp_1_time2==1)/sum(nonemp_2_time1==1);
Dni_2_Din_1=sum(nonemp_2_time1==1 & informal_2_time2==1 & informal_1_time1==1 & nonemp_1_time2==1)/sum(nonemp_2_time1==1);

transitions=[Dfn_1;Din_1;Dnf_1;Dni_1;Dfi_1;Dif_1;Dfn_2;Din_2;Dnf_2;Dni_2;Dfi_2;Dif_2;Dni_1_Dfn_2;Dni_1_Din_2;Dni_2_Dfn_1;Dni_2_Din_1];
end





